MULTIPHASE MODELING OF TUMOR GROWTH WITH MATRIX 

REMODELING AND FIBROSIS 



ANDREA TOSIN AND LUIGI PREZIOSI 

Abstract. We present a multiphase mathematical model for tumor growth which incor- 
porates the remodeling of the extracellular matrix and describes the formation of fibrotic 
tissue by tumor cells. We also detail a full qualitative analysis of the spatially homogeneous 
problem, and study the equilibria of the system in order to characterize the conditions 
under which fibrosis may occur. 



1. Introduction 

It is well known that tumor tissues are often stiffer than normal tissues. For instance, a 
normal mammary gland has an elastic modulus of about 2 hPa, which may dramatically 
increase for a breast tumor up to about 4 kPa pj. For this reason self-palpation is often 
a successful tool of pre-diagnosis for the detection of possible stiffer nodules, therefore 
so encouraged. In most cases, the increased stiffness is due to the presence of a denser 
and more fibrous stroma [21 |3l H] coming from a considerable change in the content of 
ExtraCellular Matrix (ECM). Indeed, as reported in [Ij, doubling the percent amount of 
collagen would increase the stiffness of a tissue by almost one order of magnitude (328 Pa 
and 1589 Pa for a 2 mg/ml and a 4 mg/ml collagen mixture, respectively). The percentage 
of ECM also changes within the same tumor type during tumor progression [5] . 

The continuous remodeling of ECM is a physiologically functional process, because it 
allows to keep the stroma young and reactive. In fact, prolonged rest is detrimental for 
bones and muscles, while physical training has an opposite effect. The ECM is constantly 
renewed through the concomitant production of Matrix MetalloProteinases (MMPs) and 
new ECM components. In stationary conditions, remodeling of ECM is a slow process: for 
instance, in human lungs the physiological turnover of ECM is 10 to 15% per day [6] , which 
leads to an estimated complete turnover in a period of nearly one week. However, when a 
new tissue has to be formed, e.g. to heal a wound, the rate of production is one or two order 
of magnitude faster [3 [S]. It is also well known that the remodeling process is strongly 
affected by the stress and the strain the tissue undergoes, as clearly occurs for bones, teeth, 
and muscles PdUllIlj. Hence, the relation between the rate of production/degradation of 
ECM constituents and the pressure felt by the cells is rather complicated. 
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Increased presence of ECM does not characterize only many tumors but was also ob- 
served in other pathologies like intima hyperplasia, cardiac, liver, and pulmonary fibrosis, 
asthma, colon cancer [121 HSl [HI El UHl |T6]. The alteration in the ECM composition 
can be due to several probably concurring reasons, including increased synthesis of ECM 
proteins, decreased activity of MMPs, upregulation of Tissue-specific Inhibitors of Metal- 
loProteinases (TIMPs). 

The interaction between ECM and cells is also attracting the attention of many research 
for other reasons. Indeed, on the one hand cells must adhere properly in order to survive, 
a phenomenon called anoikis, as well as be anchored to the ECM to undergo mithosis. On 
the other hand, the interaction with the stroma has been argued to be one of the causes 
of tumor progression |21 HZl |3l [151 [IB] • 

In this paper we propose a general multiphase mathematical model able to describe the 
formation of fibrosis through either excessive production of ECM or underexpression of 
MMPs. The model is based on the frameworks deduced in [191 [20] , taking also cell-ECM 
adhesion into account. In particular, ECM is regarded as a rigid scaffold while the cell 
populations (tumor and healthy cells) are assumed to behave similarly to elastic fluids. 
More realistic constitutive models, taking cell-cell adhesion into account and comparing 
theoretical and experimental results, can be found in [211 122] • For the sake of conciseness, 
we refrain from citing here all papers dealing with multiphase models of tumor growth, 
and refer to the recent reviews [231 1211 125] for more references. 

In more detail. Sect. |2] derives and describes the model, which is then studied from 
the qualitative point of view in Sect. [3| having in mind the general dependence on the 
parameters stemming from biology. Existence, uniqueness, and continuous dependence of 
the solution on the initial data is proved in the spatially homogeneous case, and equilibrium 
conflgurations are discussed. These theoretical investigations reveal several interesting 
features of the model, for instance the fact that it predicts no other equilibria but the fully 
physiological and the fully pathological ones, featuring no tumor cells and no healthy cells, 
respectively. The physiological equilibrium turns out to be stable in the manifold with no 
tumor cells, but becomes unstable as soon as few tumor cells are present, which trigger the 
formation of flbrotic tissue. 



2. Multiphase modeling: general picture and particular cases 

In the multiphase modeling approach, tumors are regarded as a mixture of several in- 
teracting components whose main state variables are the volume ratios, i.e., their percent 
amounts within the mixture. With a view to providing a simplifled, though still realistic, 
description of the system, we conflne our attention to two cell populations: tumor cells, 
with volume ratio 0t, and healthy cells, with volume ratio (pn, moving within a remodeling 
extracellular matrix with volume ratio (pM- Clearly, < 0q, < 1 for all a = T, H, M. 
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Balance equations for the cellular matter. Following [20j, we obtain the main governing 
equations for the cellular matter by joining the mass balance equation and the correspond- 
ing balance of linear momentum (with inertial effects neglected): 

where := 0r + <t>Hi is the duplication/death rate, and the cell motility tensor 
within the matrix. Cells are regarded as elastic balloons forming an isotropic fluid, and are 
assumed to feature equal mechanical properties, hence their stress tensor is T = — for 
a pressure-like function E. In addition, the model accounts for the attachment/detachment 
of the cells to/from the matrix by means of a stress threshold > 0, which switches 
cell velocity on or off according to the magnitude of the actual stress sustained by cells in 
interaction with the matrix (see [20] for further details). 

In the application to matrix remodeling and fibrosis, we consider that cells duplicate 
and die mainly on the basis of the amount of matrix present in the mixture. In general, 
also the availability of some nutrients plays a major role, but in the present context we 
assume that they are always abundantly supplied to the cells. Specifically, we set 

(2) = la{(t)M)H,^ {tpa -Sa- S'^H.^j {rria - (Pm), 

where ip '■= <Pt + (pH + 0a/ is the overall volume ratio occupied by cells and matrix, and 
7q,(-) is the net growth rate of the cell population a, tempered by the free space rate H^^. 
In particular, the H^Js are functions bounded between and 1, which vanish on (— oo, 0) 
and equal 1 on (e„, -|-oo) (further analytical details in Sect, [sj Assumption [T]) . Cell growth 
is inhibited when the amount of free space locally available is too small {ip > ipa) with 
respect to a threshold tpa £ [0, 1]. At the same time, either apoptosis or anoikis can trigger 
cell death at rates 6a, S'^ > 0, respectively, the latter taking place when a too small amount 
of ECM {(pM < ma) with respect to a given threshold nia G [0, 1] results in an insufficient 
number of possible adhesion sites. 

Usually 7t(-) = 1h{-), = ^h, = ttT'H, = = ^m- Instead, a difference between 
ipT and ipH, with tpx > ipn-, may be related to a smaller sensitivity to contact inhibition 
clues by tumor cells [19j. On the whole, we notice that it must be 

(3) Tt{(Pm, ^) > r//(0M, ^), v(0Af, ^) e [0, 1] X [0, i] 

which holds if: (i) 6'rp < 5'^ (smaller sensitivity to anoikis by tumor cells), (ii) e-r > 
(different speed for the switch mechanism, e.g. because of a different uncertainty in the 
response to mechanical stimuli), (iii) < mn (higher capability from tumor cells to 
escape anoikis by surviving a greater lack of adhesion sites). 

Matrix remodeling. In general, the ECM is a quite complicated fibrous medium. For the 
sake of simplicity, we model it as a rigid scaffold, which makes it unnecessary to detail 
its stress tensor because the internal stress is indeterminate due to the rigidity constraint. 
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Under this assumption, the evolution equation for the volume ratio 0m reads 

(4) - Tm, 

where Tm is the source/sink of ECM accounting for remodeling and degradation due to 
the motion of the cells within the scaffold. Notice that, in general, 0m depends on both 
time t and space x, although the latter acts mainly as a parameter in the above differential 
equation. 

Matrix is globally remodeled by cells and degraded by MMPs, whose concentration per 
unit volume is denoted by e = e{t, x): 

(5) Tm = ^ fia{4>M)H^j^j{i)M - i')(pa " z^e0M- 

a=T,H 

Here, /i^ is a nonnegative, nonincreasing function (cf. Sect, [sj Assumption [T]) repre- 
senting the net matrix production rate by the cell population a tempered by the free 
space rate H^j^^, and z/ > is the degradation rate by the enzymes. As usual, the 
latter are not regarded as a constituent of the mixture, but rather as massless macro- 
molecules diffusing in the extracellular fluid according to a reaction-diffusion equation: 
et = DAe + J2a=T n'^a'Pa — e/r, for net production rates tTq, > and enzyme half-life 
r > 0. Actually, enzyme dynamics is much faster than that involving cell growth and 
death, hence it is possible to work under a quasi- stationary approximation. Furthermore 
enzyme action is usually very local so that also diffusion can be neglected and fi- 
nally e = T^^^rp fjTia'Pa- Inserting this expression into Eq. ([s]) and defining Ua ■= uTTCa 
ultimately yields 



M — 

a=T,H 



^afj^M ) 00 



The pathological cases possibly leading to fibrosis are either /ir(-) > /^//(O or ut < i^h, 
which imply that tumor cells produce either more ECM or less MMPs than healthy cells, 
respectively. 

3. The spatially homogeneous problem 

The spatially homogeneous problem describes the evolution of the system under the main 
assumption of absence of spatial variation of the state variables 0t, 4>h, 4>m- This allows 
in particular to describe the equilibria, and the related basins of attraction, as functions 
of the parameters of the model. 

In the sequel we will be concerned with the following initial value problem: 



(6) 



dt 



Yfc{<pM)H^^ {tpa -ij) -Sa- 5'o,HeM i^a " 0J\/)] 0q, a = T, H 



dt ^ 

a=T,H 

«(0) = 0° G[0, 1], a = T,H,M 
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over a time interval (0, T], T > 0. Some preliminary technical assumptions are in order: 

ASSUMPTION 1. For a = T, H, M as appropriate, we assume < ipa, f^a < 1; 
6a, S'a, Va > 0, and in addition that 7q,, /Iq, : [0, 1] IR+ are Lipschitz continuous, with 7^, 
nondecreasing, 7^(0) = 0, and Ha nonincreasing. 

Moreover, we assume that the functions H^^ : M [0, 1] are Lipschitz continuous and 
vanishing on (—00, 0]. 

The monotonicity of 7q,, is dictated by the fact that cell proliferation is fostered by 
the presence of ECM, whereas production of new matrix is progressively inhibited by the 
accumulation of other matrix. Owing to the properties recalled in Assumption [1} 7q,, /Xq,, 
and H^^ satisfy 

(7) 7a(s) < Lip(7a)s, 7„(s) < 7a(l), < < /ia(0), Vs G [O, l] , 

(8) <Lip(i7,J|s|, VseM, /3>0. 

The functions H^^ may be taken to be mollifications of the Heaviside function, for instance 
= if s < 0, H^^{s) = e~^s if < s < and H^^{s) = 1 if s > or even 
smoother. 

Let us introduce the space := C([0, T]; M'^) of continuous functions u : [0, T] M'', 
endowed with the 00-norm ||u||oo = maxtg[o,r] In proving our results we will utilize 

d = 3 and 4 = 4. 

Well-posedness. We start by studying existence, uniqueness, and continuous dependence 
on the data of the solution = (0^, 4>h, <Pm) to problem (|6|. We will then also discuss 
its regularity. Since the 0q's are volume ratios, we are interested in nonnegative solutions 
such that ^(t) < 1 for alH > 0. 

THEOREM 1 (Existence, uniqueness, and continuous dependence). For each initial da- 
tum 4>^ > with < 1, problem (|6| admits a unique nonnegative global solution 
4> € C([0, +00); M^) such that ||0||oo < 1- I'l^ addition, if <pi, 4>2 are the solutions corre- 
sponding to initial data 0°, cf)^, then for each T > there exists a constant C = C{T) > 
such that 

||02-0llloo<C(T)||0°-0?||i 

in the interval [0, T]. 

Proof. 1. Let us introduce the function ip{t) := 1 — ip{t) (which, in mixture theory, 
identifies the free space available to be filled by some extracellular fluid) and consider 
the auxiliary problem given by the set of equations ^ plus = — J2a ^'a^ along with 
(f^ := (p{0) = 1 — J2a't'a- Clearly, a triple (0t, 0Af) is a solution to problem ^ if and 
only if the quadruple (0t, (pH, (pM, f) is a. solution to the auxiliary problem. 

2. We put the auxiliary problem in compact form as 

— = J * , t > 

(9) I dt ^ ^' 

I $(0) = 
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where $ = (0, if) and J : V"^ —>■ is given componentwise by the right-hand sides of the 
differential equations in ^ plus = — J2a '^a- 

Next we make the substitution = ^{t)e~^^, where A > will be properly selected. 
Due to the specific expression of J, the term J[^(t)e~^*] can be given the form /[^'](t)e~'^* 
for a suitable operator I : ^ , which allows us to rewrite problem ^ in terms of ^ 
as 

^ = /[*] + A*, t>0 
*(0) = $° 

or, in mild form, as = $° + / (l[^]{r) + A*(r)) dr =: G[^]{t). 

^ ^ 

3. Let us look for a mild solution in the following set of admissible functions: 

A = {ueV^ : u{t) > 0, ||u(t)||i = e^* for all t G [0, T]}. 

Notice that & A amounts in particular to (pa(t), ip(t) > with J2a 'Paif) + ^if) = 1 for 
all t G [0, T], thus < 1, which is what the saturation constraint requires on the 

volume ratios of the constituents of the mixture. 

Any mild solution of ^(t) = is a fixed point of the operator G, therefore the 

task is to show that G admits a unique fixed point in A. 

4. Owing to Assumption [l] and properties ([7|, ([s]), if u(t) > then 

t 

Go\n]{t) > 0° + (A - C„) j u^{t) dr {Ct^h = 5t,h + 5t,h, Cm = i^t + t^h), 



G'^[u](t)>^°+ (A- (7a(l)Lip(i^eJ+/ia(0)Lip(i/,,J)) /"«^(r)rfr, 

\ a=T,H / { 

hence we can choose A > so large that G'[u](t) > as well. If in addition ||u(t)||i = e^^ 
then, using = — ^„-^a, we discover ||G'[u](t)||i = e'^*. In conclusion, u E A implies 
G[u] G A, i.e., G maps A into itself. 

5. Take now u, v G ^ and observe that 

t 

|i + A||u(r)-v(r)||i) dr. 



||G[u](t)-G[v](t)|K< J (||/[u](r)-/[v](r) 







Lipschitz continuity of 7^, jja, H^^ along with H^^(s) < 1 and properties ([T]), ([s]) imply 
that there exists G > 0, independent of T, such that |Ja[u](t) — /Q[v](t)| < C||u(t) — v(i(:)||i 
each a = T, H, M. Since = — Xla-^"' analogous relationship holds true also for 
I^, hence finally \\G[u] — G[v]||oo < T{G + A)||u — v||oo, which proves that G is Lipschitz 
continuous on A. 
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6. From the above calculations we see that we can choose T > so small that G 
be a contraction on A. Since ^ is a closed subset of V^, Banach Fixed Point Theorem 
asserts that G has a unique fixed point ^ & A. Therefore, the auxiliary problem ^ 
admits a unique nonnegative local solution $ G V^'^ such that ||$(t)||i = 1. The three 
first components of $ form the unique nonnegative solution (j) eV^ to problem ^ with 
ll0WI|i<l. 

Next, taking 0(T) as new initial condition and observing that it matches all the hypothe- 
ses satisfied by 0°, we uniquely prolong over the time interval [T, 2T] in such a way that 
4>{t) > and II (/)()!:) Ill < 1 all t G [0, 2T]. Proceeding inductively, we ultimately end up 
with a unique nonnegative global solution G C([0, +oo); M^), for which the estimate 
||(^||oo < 1 easily follows from ||0(t)||i < 1 all t > 0. 

7. Let now ^'i, ^2 G ^ be the two mild solutions corresponding to initial data $2- 
Using the previous estimates we discover 

t 

||*2(t)-*iWI|i < ||$^-$?||i + (C + A) j ||*2(r)-*i(r)||irfr, 



whence, invoking Gronwall's inequality, 

||*2(t) - *i(t)||i < [1 + (C + A)te(^+^)*] ||$° - $?||i. 

Returning to $1, $2 and observing that ||$2 ~ ^illi ^ 21102 — ^'j'Hi we finally get the 
desired estimate of continuous dependence, after taking the maximum of both sides for 
tG[0, T]. □ 

THEOREM 2 (Regularity). // the functions 7q,, /Iq,, H^^ are of class on [0, 1] then 
the solution (p is of class C^^^ on [0, +00). 

Proof. According to Theorem [l| the 0a's are continuous on [0, +00), therefore the right- 
hand sides of the differential equations in ^ define continuous functions on [0, +00). It 
follows that the (p'^^s are continuous as well, i.e., the solution is actually on [0, +00). If 
7a, /ia, H^^ are of class C'' then, by differentiating the ODEs in (|6| /c times, this reasoning 
can be applied inductively to discover G C'^+"^([0, +00); M^). □ 

Stability of the equilibrium configurations. Next we study the equilibria of model (|6|. It 
is immediately seen that (pT = 4'h = ^ gives rise to an equilibrium solution for any 
(pM £ [0, 1], corresponding to that all cells have died leaving some ECM. In order to 
investigate nontrivial equilibrium configurations, we proceed by considering first the two 
important sub-cases in which either 0t = or (pn = but (pT, (pH do not vanish at the 
same time. The former will be referred to as the fully physiological case, the latter as the 
fully pathological one. In the following, (pa will be the nonzero volume ratio for either 
a = T 01 a = H , meaning that the other one is identically zero. For the sake of simplicity, 
let us fix %pM = 1 and examine the case (pa + (pM ^ — V for some arbitrarily small ?7 > 0, 
in such a way that, choosing < V: "we have the simplification H^^j{%pM — (pa — (pu) = 1- 
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Figure 1. Left: phase portrait in the fully physiological or pathological 
case. Right: phase portrait for the full model. 



(10) 0A/ = M„, 0, = - - h: 



Suppose that the function fia{s) — VaS has exactly one zero, say s = Ma- Since, for 
physiological reasons, we further must have fiai'ipa) < ^ai^a, it follows Ma G (0, ipa)- In 
this case, there is one nontrivial equilibrium given by 

7a (M„) 

which is readily checked to be stable. Notice that the function H~^{s) is well defined for 
s E (0, 1). 

The trivial equilibrium with also 0q = can be reached basically in two situations. The 
first one is when (pM is initially too small, so that the growth rate of the cells is lower than 
the apoptotic rate and anoikis occurs. This corresponds to initial conditions located in the 
lower-left corner of the phase portrait illustrated in Fig. [l| left. The equation = (/>a((^M) 
of the curve delimiting this basin of attraction in the phase space is obtained by integrating 

^^^^ d(j)a _ la{4>M)H^^{llja - (pg - 0m) - Sg - dgH^^jimg - (Pm) 

d(f)M IJ'a{(f)M) - l^a(t)M 

with the condition (pg{4>*Mg) = 0, (p^.^^ ^ (O5 1) being the smaller root of the equation 
Fa (0m) = with (pT = (f>H = 0, cf. Eq. ([2]), characterized by F'„(0^j^) > 0. This region 
does not exist if (pXjg < 0. 

The second case is when (pM is initially too large, namely the ECM is overly dense and 
cells are so compressed that the growth rate is again lower than the apoptotic rate because 
Heaii^a — 4'a — 4'm) ~ 0. This corrcspouds to initial conditions located in the lower-right 
corner of the phase portrait depicted in Fig. [T| le ft. The curve delimiting the basin of 
attraction is again obtained by integrating Eq. now with the condition 0«(0mq) = 0, 
4>Ma ^ (O5 1) being the larger root of the equation Tg{(f)M) = with 4>t = 4>h = 0, so that 
0Mo — 4>Ma- this case F'„(0^„) < 0. The region does not exist if > 1. 

In order to get the complete picture, we further have to investigate whether a nontrivial 
equilibrium with 0^, (pn > may exist. For this, we recall that the duplication/death rates 
Tg are constructed so as to match the biological requirement Tt{(Pm, i') > ^nifpAi, "0) for 
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Figure 2. Left: formation of normal tissue in the physiological case. Right: 
formation of hyperplastic fibrotic tissue due to a small initial amount of 
tumor cells. 



all (0M, i^) e [0, 1] X [0, 1], cf. Eq. (§. As a consequence, we see that it is impossible for 
the right-hand sides of the two first equations of problem ^ to vanish simultaneously at 
an equilibrium point (0t, 4>h, 4>m) with 0^, (pn 7^ 0, for this would imply that there exist 
4>M £ [0, 1] and ip ^ (0, 1] such that Ta{4>M, V') = for both a = T, H, which contradicts 
the above-mentioned Eq. ([3]). 

Hence, the only possible equilibria of the system are those arising in the fully physiolog- 
ical or pathological situation. In addition, condition ([s]) makes the nontrivial physiological 
equilibrium unstable and the nontrivial pathological one stable, as it can be realized from 
the three-dimensional phase portrait shown in Fig. [ij right. 



Figure |2] (left) shows an example of a temporal evolution of the system giving rise to the 
formation of normal tissue in the fully physiological case. The initial death of healthy cells 
is due to anoikis, indeed cells are seeded in an environment completely deprived of ECM, 
that they have to build fast enough. The decrease stops as soon as the amount of ECM 
produced is such that 'yH{<pM)H^^{ipH — ip) > Sh + ^'H^eMi^H — (pAi), then the number of 



cells starts increasing, eventually leading to the stationary solution predicted by Eq. (10) 
for a = H. Conversely, if the initial amount of cells is insufficient to produce ECM rapidly 
enough then the entire population will die. 

Figure [2] (right) gives instead an example of a complete temporal history ending with the 
formation of hyperplastic and fibrotic tissue. Despite the initial conditions 0^, coincide 
with the equilibrium values reached after the formation of normal tissue, the presence of a 
small amount of tumor cells {(f)^ > 0) changes dramatically the outcome of the evolution, 
leading to a full depletion of healthy cells. This evolution can be duly compared with that 
shown in Fig. |3} which refers to the simulation of the full spatial and temporal model in 
one space dimension, cf. Eqs. ([T]), (|4]). Starting from the same initial conditions as in the 
spatially homogeneous case, the presence of a small amount of tumor cells at the beginning 
generates a traveling wave, which progressively depletes healthy cells and produces further 
fibrotic matrix while invading the normal tissue. 
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Figure 3. Traveling wave solutions for <^ii (blue), 0t (red), 0m (green) in 
the full spatial and temporal evolution. 
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